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We derive the density functional for the ground-state energy of a two-dimensional, spin-polarized 
gas of neutral fermionic atoms with magnetic-dipole interaction, in the Thomas-Fermi-Dirac ap- 
proximation. For many atoms in a harmonic trap, we give analytical solutions for the single-particle 
spatial density and the ground-state energy, in dependence on the interaction strength, and we dis- 
cuss the weak-interaction limit that is relevant for experiments. We then lift the restriction of full 
spin polarization and account for a time-independent inhomogeneous external magnetic field. The 
field strength necessary to ensure full spin polarization is derived. 
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I. INTRODUCTION 

U ltracold atomic gases provide highly controllable sys- 
tems for the study of condensed-matter phenomena 
With the ongoing experimental efforts in ultracold gases 
of neutral Fermi atoms [2J and the possibility of a gen- 
uinely two-dimensional (2D) confinement [3t3I, it is now 
possible to acquire data about dilute 2D degenerate 
Fermi gases of neutral atoms under adjustable labora- 
tory conditions [lH • It is hoped that this will advance our 
understanding of various 2D phenomena, such as high-T c 
superconduction || Q , effective massless Dirac fermions 
[13, [Hi, the BEC-BKT cross over [HE3, and others - 

While experiments of this kind will involve periodic 
potentials of various geometries, all experiments with ul- 
tracold atoms in a 2D configuration will make use of a 
confining potential. Before investigating aspects of spe- 
cific periodic potentials, one can examine the behavior 
of the degenerate gas in the 2D trapping potential. Wc 
study this situation with the help of density function- 
als, and arrive at detailed predictions for the case of an 
isotropic harmonic trapping potential. 

Density functional theory (DFT), which has its histor- 
ical roots in the Thomas- Fermi model for atoms [l5|, [l6[ , 
was first formulated for the inhomogeneous electron gas 
[T3|, with immediate applications to atoms, molecules, 
and solids; see Ref. 18], for example. DFT can equally 
well be used for studying other physical systems, such as 
dilute gases of neutral fermionic atoms under the influ- 
ence of a confining external potential. While the DFT 
formalism can be based on both the spatial [13] and the 
momental density (l9| . the spatial-density version gives 
a more natural description in the case of a position- 
dependent interaction, such as the magnetic dipole in- 
teraction. We derive the density functionals for spin- 
polarized fermions with magnetic dipole interaction, con- 
fined in a 2D harmonic potential, and investigate the 
ground-state density and energy of the system. 

The article is organized as follows. Section [TT] sum- 
marizes earlier investigations in three dimensions (3D). 
In Sec. IIII1 we discuss how to appropriately reduce the 



dimensionality. Various 2D density functionals are de- 
rived in Sec. lIVI A discussion about the scaling behavior 
of these functionals is then given in Sec. [V] Section IVII 
presents the analytic results of the ground-state density 
and energy, and discusses the weak interaction limit. Sec- 
tion [VlT] extends the formalism to accommodate the spin 
dependence in an inhomogeneous magnetic field. We con- 
clude with a summary and a brief outline of prospective 
work. 



II. THE 3D CASE 

It is expedient to recall some basic relations that were 
earlier established in 3D, mainly collected from Refs. [20l — 

m. 

A. Single-particle density and density matrix 

The spatial one-particle density matrix n^(r'; r") and 
the associated one-particle Wigner function z/(r, p) are 
related by 

(1) 

with (dp) = dpxdpydp z denoting the volume element in 
the momentum space. The spatial and momental one- 
particle densities are obtained by integrating u(r, p) over 
the other variable, 

„(r) s „<«(r ; r) = JMj^^p), 

'<"> = / " (r ' P) ' <2> 

Note that both densities are normalized to the total num- 
ber of particles N, 

N = [ (dr) n(r) - / (dp) p(p) . (3) 
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B. Density functionals for energy 

For a system of spin-polarized fermions in an isotropic 
harmonic trap, the potential energy is given by 



£tra P N = I (dr)^Muj 2 r 2 n(r), 



(4) 



where M and w are the mass of an individual atom and 
the trap frequency, respectively, r = |r| is the length of 
the position vector r, and the kinetic energy is 



Skin = J (dp) iL p( P ) . 



(5) 



Both Strap and Skin are sums over single- particle contri- 
butions. 

The interaction energy, E dd , which is a sum 
over particle-pair contributions, is evaluated using 
the diagonal part of the two-particle density matrix 
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E dd = - I (dr')(dr")V dd (r' -v")n^(r\r";r',r") (6) 



with the magnetic dipolar interaction potential 

,2 /,. „\2 



(7) 



where fi and /x are the magnetic dipole moment and its 
magnitude of an individual atom. The contact term in 
Vdd (r) is necessary to ensure that the magnetic field gen- 
erated by the point dipole is divergence-free. 



one has to use v(r, p) = r)(t(p) — V(r)) where V(r) is 
the external potential and t(p) is determined by p(p) 
through Eq. ([2]). In the case of an isotropic harmonic 
potential, V(r) cx r 2 , this is 



<(r,p)=T 1 (h[6n 2 p(p] 



1 1/3 



(10) 



visibly the analog of Eq. ©. 

This yields the familiar density functional of the kinetic 
energy, 



SkinW 



ft 2 1 

(dr) M 20n< 



[67T 2 n(r)] 



5/3 



(11) 



Since the contributions of the contact term to the direct 
and to the exchange energy cancel each other in the fully 
spin-polarized situation under consideration, and the re- 
maining exchange energy vanishes under the average over 
the solid angle associated with the relative distance, the 
density functional of the interaction energy, 



E dd [n] 



1 



(dr)(dr>)n(r)V dd (r - r>)n(r>) , (12) 



is characterized by an effective potential V dd , 

>-r) 2 " 



(13) 



D. Ground-state energy and density 



C. TFD approximation 

In the spirit of the approach that was pioneered by 
Thomas, Fermi, and Dirac (TFD), a two- fold semiclas- 
sical approximation is employed here. First, nS 2 ' is re- 
placed by products of factors (due to Dirac (24|) 
according to 

n(Vi,4;r?,r 2 ') = n^(r[;r'{)n^(r' 2 ;4) 

-n(Vi;r> ( V 2 '; I"'/)- (8) 

This splitting corresponds to the direct and exchange 
terms when evaluating the interaction energy, E dd . Note 
that this expression is only valid if the system is spin- 
polarized. Otherwise, a multiplicative constant preced- 
ing the second term is needed to account for the spin- 
multiplicity. 

Second, the Wigner function is a uniform sphere of a 
finite size (due to Thomas [l5[ and Fermi [l6j ) 

i/(r,p) = V (h[67r 2 n(v)] 1 / 3 -p), (9) 

where rj( ) is the Heaviside unit step function. This ap- 
plies when functionals of the spatial density n(r) are con- 
sidered. For functionals of the momental density p(p), 



Thus, the functional for the total energy of the ground 
state in the TFD approximation is given by the sum of 
the three terms in Eqs. g]), (HU), and (fl"2]) . 

S TPD [n] = |(dr)|^L[ 67 r 2 n(r)]^ 
+ J(dr)^Muj 2 r 2 n(v) 

+i J{dr){dr')n{r)V Ad {v-r')n{r'). (14) 

Upon applying the variational principle, we find that 
the density that minimizes -Etfd must obey the integral 
equation 

2M 1 WJ 2 

+ j{dv')V dd {v-r')n{v') = X -Muj 2 R 2 , (15) 

where ^Moj 2 R 2 is a convenient way of parameterizing 
the Lagrange multiplier for the normalization constraint 
ofEq. ©. 
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III. FROM 3D TO 2D 



C. Various energy terms 



The form of the density functional in Eq. (114)) gives 
no explicit indication of its dependence on the spatial 
dimension. It is thus necessary to re-derive the density 
functionals in 2D, with some suitable assumptions about 
the Wigner function. 



A. A possible Wigner function 

When the trapping potential in the z-direction is har- 
monic and sufficiently stiff, as is the typical situation 
in an actual experiment, the system will remain in the 
ground state in this direction, and this gives rise to a fac- 
torizable Gaussian dependence in z and p z in the Wigner 
function, 



pi = i'_(r± . I' i 2 (-xi) j ~ - ^ 



06) 



where l z = ^Jh/(Muj z ) ' is the harmonic oscillator length 
scale in the z-direction, and the subscript indicates 
that these various quantities live in the transverse xy- 
plane. Here, ui z is a finite but large frequency and, in or- 
der to achieve a 2D geometry, we require that fkd z ^> ksT 
for the situation of ultracold atoms that we have in mind, 
although we take the limit uj z — > oo for mathematical 
convenience whenever possible. 



B. Densities in 2D 

In analogy with the densities defined in 3D, Eqs. ^ 
and © , the densities in 2D are given by, 



{2irh) 2 

, x f (dp±) i , 
n ^ = J (2W Mr± ' P±) ' 
f (dr±_) 



(17) 



With the decomposition of the Wigner function in 
Eq. ([To]) , we find that the densities in 2D and those in 
3D are related in the following manner: 



n<V;r") =^ 1) (rl;r'|)— ^ exp 
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n(r)=^(r ± )-^exp(-|), 



p(p) = PX(P±) ^ exp , (18) 

where z + — i(z' + z"), and z_ = z' — z", such that the 2D 
densities are now normalized to the number of particles, 

N = J(dr 1 _)n x (r 1 _)= f(dp±) p±(px) ■ (19) 



By integrating over z and p z , we immediately find the 
trap energy and the kinetic energy in terms of the 2D 
densities, 

Strap = J (dr x ) nx(rx) \^Mu\r\ + ^hu z , 
f v 2 N 

S ki „ = / (dp±) P±(P±) ^ + ~^z , (20) 

where lu± is the radial trap frequency in the xy-plane, 
assuming isotropy. Note that both expressions contain 
parts analogous to the corresponding expressions in 3D 
and additive constants, which are the sum of single- 
particle energies in the ground state of the harmonic trap 
of the z-confinement. Since these constants play no role 
in the dynamics of the system, we renormalize the ex- 
pressions, such that 



Strap — 



(dr ± )n ± (r x )^Mulr 2 ± , 



f P 2 
Skin = (dp±)p±{p±) ^ , 



(21) 



which are now independent of lo z and unaffected when 
the limit u z — > oo is taken. 

To investigate the interaction energy, we employ the 
(2+l)D version of Eq. ©, 



n( 2 )(r',r";r',r") = n(r')n(r")-n (1) (r';r")n (1) (r";r') 



r") - nF> 



l 2 ir 



(ni(^)ni(rl) 

-n^^yDn^Ksri)), (22) 

which corresponds to the splitting of the direct and ex- 
change energies. It is clear from Eq. (|22|) that the con- 
tact term in the interaction potential, Eq. (J7J), enforces 
r' = r" and thus equal and opposite contributions from 
the direct and exchange energies, the familiar situation 
when the system is spin-polarized. It is then permissible 
to drop the contact term, which amounts to replacing the 
original interaction potential by the effective potential of 
Eq. (USD, i.e. 



S 



dd 



(dv')(dr")V dd (r' -v") 
fn(r')fi(r")-n (1) (r';r> (1) (r"|r')) • (23) 



Since V dd (r' — r") depends only on the difference in the 
positions, we isolate the z-direction by identifying p = 
(r'-r"U, 



V dd (r'-r") = 



47T 



3z_ 



{p 2 + z 2 f/ 2 (p 2 + z 2 fl 2 



(24) 
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where we assume that the magnetic dipole moments of 
all fermions are polarized in the z-direction, i.e. fi = fj,e z . 

It should be noted that the replacement of Vdd( r ) 
by V^dd(r) only takes place after the approximation in 
Eq. (|22p is made. It may not be correct for a differ- 
ent approximation scheme, that is: when going beyond 
Dirac's approximation in Eq. ([5]). 

In the limit of u> z —> oo, the Gaussians of z± become 
Dirac delta functions, yielding 



E, 



del 



- (dr ± )(dr ± ) — 



1 



x (niK)nL (r'i ) - n? } (r' ± ; r'{ )n« (r'| ; r' ± ) 

(25) 

In hindsight, we recognize the result above as an imme- 
diate consequence of having fi _L r, which forces their 
scalar product in Eq. (fT3| to vanish, while the contact 
term does not contribute for the reason discussed earlier. 



IV. 2D FUNCTION ALS 

From this section onwards, for notational convenience, 
we leave out all the subscripts It is understood that 
all the densities refer to the 2D definition specified in 
Eqs. (fT"7|) . and all the vectors reside in the xy plane. 

A. TFD: A brutally simple Wigner function 



3D functional in Eq. (TTT1) . and the ID functional for the 
kinetic energy has the cube of the density; see Table U 
below. 

The interaction energy, in particular, turns out to con- 
sist of two pieces with different dependence on the one- 
particle density, namely 



^dd - ^dd 



E 



(2) 
dd > 



47T 



(dr)^V^n(r) 5 / 2 . 
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#ddM = -^-tt /WnWv^nfr) , (29) 



where \J — V 2 ' is an integral operator that is given by 



V /= V^n(r) = J^(dk)k 



n(r') . (30) 



We report the details of deriving Eqs. (f2T)f in the Ap- 
pendix. 

Note that the splitting of the interaction energy in 
Eqs. (|2"9")l does not correspond to the direct and exchange 
energies as in the 3D case, where, as we recall, the ex- 
change energy exactly compensates for the contribution 
of the contact term to the direct energy. In 2D, both the 
direct and exchange energies are infinite individually, and 
they can only be considered together so that the total in- 
teraction energy is finite. Both contributions in Eqs. (|29[) 
stem from the sum of the direct and the exchange energy. 



In order to derive the density functionals in 2D, one 
first needs to complete the TFD approximation started 
in the previous section and assume a "brutally simple 
Wigner function" '23| analogous to Eq. ©, i.e. 



f(r,p) = ?y(/i[47rn(r)] 1//2 - p) 



(26) 



where the power and prefactor of the density are deter- 
mined by normalization. 



B. Density functionals 

Not surprisingly, the potential energy reads 



£tra P M = J (dr) ilfwV n(r) , 



(27) 



where we emphasize that uj is the radial trap frequency in 
the xy- plane, assuming isotropy. Upon using Eqs. (|I7p. 
(|2ll. and (J26J) , we find that 



E kin [n]= J(dv)^nn(r) 2 



(28) 



The n(r) 2 dependence of this 2D functional can also 
be obtained from dimensional analysis; similarly, dimen- 
sional analysis confirms the n(r) 5 / 3 dependence of the 



C. Ground-state energy and density 

The total energy of the system in the TFD approxima- 
tion is now given by the sum of the various energy terms 
derived above, 



E TFU [n\ 



47T 



h 2 r l 

(dr) — 7rn(r) 2 + (dr) -Muj 2 r 2 n(r) 
f(dr) —^Fn(rf' 2 - 7rn(r) v / -V 2 'n(r) 

J L 45 



(31) 

The density that minimizes the energy, constrained by 
the normalization condition (|I9j) . must then obey 



2h , N 1 » r 2 2 

7m(r) + -MwV 

M w 2 
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4tt 

27iV-V 2 ' n(r) 



^^n(rf' 2 



-Muj 2 R 2 . 



(32) 



where, as in Eq. (TT5l . ^Muj 2 R 2 is the chemical potential. 

A comparison between Eqs. (f3"2"|) and (Tl~5|) shows that 
the reduction of dimension does not yet provide any op- 
erational simplification when it comes to solving for the 
spacial density, because of the occurrence of the integral 
operator \J — V 2 '. But one should not fail to notice that 
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the interaction is now made up of two contributions with 
different dependence on the density. As we will see in 
the next section, the integral term is rather unimportant 
in certain parameter regimes of interest and can then be 
neglected. 



V. VIRIAL THEOREM AND SCALING 
A. Scaling transformation 

Let us consider scaling transformations that change 
both the length scale and the number of particles, 

n(r) X 2+a n(Xr) , N -> X a N . (33) 

They are consistent with the normalization constraint, 
Eq. (Ti~9")) . and affect the various terms of -EpFD in the 
following manner: 



^kin 


. \ 2+2a 17 
— > A -C/kin ) 




-^trap 


— > X 2+a Et rap , 




E {1) 


A "^dd ' 




^dd 


A Aid > 


(34) 



so that the total energy E = -E4td transforms in accor- 
dance with 

E = -Ekin + E tTa ,p + E^J + E^ 

— > A 2+2Q _Ekin + A 2+ " Strap 

+ X 3+5a ^ 2 E^J + X 3+2a E { ^ . (35) 



B. Virial theorem 

Since the minimum of E is achieved by the true 
ground-state density, all first-order changes of E in the 
vicinity of A = 1 must be generated by the explicit change 
in N, 5N = SX aN, so that 

QE 

aN dN = (2 + 2a)jBkin + + a ) E ^p 

+ (3+|«)4 1 d ) + (3 + 2a)4 2 d ) (36) 

is true for all values of a. Choosing two values of a for 
independent statements, we have 

2E kin - 2£ trap + 3 (e£> + E (2 J ) = (37) 
for a — 0, and 



on [i, u), and M 
d_ 

d 



^E = 2E dd = 2(E^+E$), 



M-Q-j^E = Strap — Skin • 



(39) 



Now, owing to the scaling argument, which will be pre- 
^dd/^dd ~ which allows 



sented next, we find that E^J/E^ 
us to neglect E^j for large A. Applying E^ ± E^J « 

E dd t0 the nrst line of E( l- <ESJ) and Ec l s - (EH)-© then 
yields 

E(jjl, u), M, A) w fuvN 3/2 £ (eA 1/4 ) , (40) 
'(M (41) 



where 



47rZj] . 



is a dimensionless interaction strength that can be un- 
derstood as the ratio between the interaction energy of 
two magnetic dipoles separated by Iq = y/h/(Mid)' and 
the transverse harmonic oscillator energy scale, and £ ( ) 
is a dimensionless function of eA 1 / 4 . We remark that 
the A 3 / 2 dependence in the prefactor results from the 
degeneracy of the harmonic confinement in 2D. 



C. Dimensionless variables 

We define the natural length scale of the system, a, the 
dimensionless position, x, and the dimensionless density, 
g(x), in accordance with 

a = l N 1 /\ x=-, 5 (x) = ^n(r), (42) 

so that the scaled density is normalized to unity. Choos- 
ing hjjN 3 / 2 as the energy unit, we have 



^tfd [g] 
fuvN 3 / 2 



.dE 



2E kin + fOStrap + (4d - Ai d ) = 4A— (38) 
for a = — |, which are supplemented by the first line of 

Eq. (|3"5")l . Further, we note the parametric dependence glected for large A on the way to Eq. (|4T)]) . 



tt |(dx) 5 (x) 2 + i J{dx)x 2 g{-x) 

+ £ iV 1 /4^^| (dx)5(x) 5/2 

- A^tf y (rfx) .g(x)V-V 2 '.g(x)^ , (43) 

where —V 2 now differentiates with respect to position x, 
and the scaled density that minimizes E^Q must obey 

2.g(x) + ix 2 + £ A/ 1 /4^ VFff(x) 3/2 

- A- 1 / 2 2 7 rV-V 2 '.g(x)^ = ^A 2 , (44) 

where \X 2 is the scaled chemical potential. The term 
preceded by A -1 / 2 originates in E?), which was ne 
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1 8 1/4 2 
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FIG. 1: The dimensionless spatial density g(x) at various 
values of eiV 1/4 = 0.01,0.1,1,10 (thin lines). The TF pro- 
file (thick dashed line) is included as a reference. Note that 
there is an insignificant difference from the TF profile for 
dV 1 / 4 < 10~ 2 . 
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FIG. 2: The energy contributions (a) £ki n , (b) Strap, and (c) 
(in units of fkuN 3 ^ 2 ) as functions of eN 1 ^ 4 (in logarithmic 
scale). The solid lines show the full solutions of Eqs. (|48[) and 
(|49p . and exhibit the correct weak- interaction of Eqs. (|50[) on 
the far left. The short-dashed lines indicate the asymptotic 
forms of Eqs. (|51[) for large values of eiY 1 / 4 . 



VI. DENSITY AND ENERGY OF THE 
GROUND STATE 



For N ~ 10 4 , which is a modest value for typical ex- 
periments with ultracold atoms, the TV -1 / 2 term is a cor- 
rection in the one-percent regime. Given that the TFD 
approximation is generally introducing errors of the order 
of a few percent, this term is of a negligible size. There- 
fore, we shall consistently discard it and all other N^ 1 / 2 
terms. Equation (131)) then reduces to 
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1 



where R TF = V2'{2N) 1 /% is the Thomas-Fermi radius 
in 2D. 

To evaluate the ground-state energy, we recognize that 
Eq. ((45]) provides a natural way of changing the integra- 
tion variable, 

-xdx=(Kg 1 / 2 + 2n)dg with k = eN 1 / 4 ^-^ , (47) 
where the position dependence of g is left implicit, such 

(2) 

(X 2 —x 2 ) (45) that au terms except E^ can be analytically expressed 



which does not single out any spatial direction and thus 
implies an isotropic ground state density, <?(x) = g(x). 
We also recognize that Eq. (j45)) is a cubic equation for 
\J g{x)' and can be solved analytically. 

In Fig. [1] we plot the dimensionless density g(x) for 
different values of eTV 1 / 4 . We observe that the stronger 
the dipole repulsion (larger e) , the lower the central den- 
sity and the larger the radius of the cloud. This feature 
is reminiscent of the condensate wave function of bosonic 
atoms when a repulsive contact interaction is taken into 
account in the mean- field formalism 25]. In contrast to 
that exhibited by a 3D spin-polarized Bose-Einstein con- 
densate [111, the simple symmetry of the isotropic har- 
monic confinement is preserved in the ground-state den- 
sity in 2D. We remark that this is partially a consequence 
of choosing the direction of spin polarization along the 
z-axis. The situation is markedly different, and more 
interesting, when the polarization direction breaks the 
axial symmetry. This will be discussed in Sec. IVIII 

On the other hand, the dipole interaction for alkali 
metals are typically small. In the limit of e — > + , we re- 
cover the well-known Thomas-Fermi (TF) profile of non- 
interacting fermions in a 2D harmonic trap, 



n(r) 



TF 



r 2 ) for < r < R TF , (46) 



in terms of k and the central density G = g(x = 0), 
Ekin = l^( 3K G 7 / 2 + 77rG 3 



hwN 3 / 2 21 



trap 



- = \x 2 - f (7k 2 G 4 + 40™G 7 / 2 + 56tt 2 G 3 ) 



huiN 3 / 2 2 21 

= J^( 7k 2 G 4 + 16tt K G 7 / 2 ) . (48) 
HujN 3 / 2 105 v ; ' 

The values of G and X are in turn determined by 



kG 3 ' 2 + 3ttG = -X 2 , 
4 

2kG 5/2 + 5ttG 2 = — , 

27T 



(49) 



of which the top equation is Eq. (|45)) for x = 0, and the 
bottom equation states the normalization of <?(x) to unit 

integral. The analytic solutions for Ey_ ini E tl:ap , and E^J 
as functions of eiV 1 / 4 = 3k/ (64^/^) are plotted in Fig. [21 
On the far left in the figure, we observe the equipartition 
of kinetic and trap energies at vanishing interaction as 
one expects in the case of a harmonic trapping potential. 

For weakly interacting atoms, i.e. eiV 1 / 4 <C 1, we ob- 
tain the various contributions to the energy up to the 
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first order in eiV 1 / 4 



fiwiV 3 / 2 3 105tt 

Strap = Vj[ 128 pl/4 yl/4 
frwiV 3 / 2 3 1057T 



1) 
dd 



_ i^-2 1 /*dV 1 /*« 0.615 dV 1 /*. (50) 
huN 3 / 2 315?r v y 

Note that the sum of Sdn and Strap has no first-order 
correction. 

The asymptotic values in the limit of large eN 1 / 4 - 
shown as dashed lines in Fig. [5] — are given by 

Skin 



TP 



hujN 3 / 2 ' tnoN 3 / 2 



S 



(2) 
dd 



( £iV l/4 ) -2/5 ^ 
(^1/4)2/5 ; 

(eJV 1 / 4 )" 1 - 58 , 



fe„W~ , (51) 

where the final power law is obtained by a numerical 
fit. Note that S trap and have the same large-eA 1 / 4 
behavior. 



VII. SPIN-DENSITY MATRIX 

While the above treatment yields the TFD approxi- 
mated ground-state density profile and energy for a 2D 
cloud of spin- 1/2 fermions that are polarized along the 
axial direction and are hence repelling each other, the 
lack of spherical symmetry of the magnetic-dipole inter- 
action, which is the source of some interesting predictions 
26] , is not well reflected due to the peculiarity of both 
the configuration and the low dimension. 

In order to take the spin-dependent nature of the 
magnetic-dipole interaction into consideration, we ex- 
tend the treatment by (i) introducing an external mag- 
netic field strong enough to define a local quantization 
axis; and (ii) constructing spin-dependent Wigner func- 
tions and hence the corresponding one-body and two- 
body spin-density matrices. 

For an arbitrary time-independent external magnetic 
field, 

B(r) = B(r) e(r) , (52) 
the magnetic energy of a single dipole is given by 

- B(r) • n = -B(r)(i e(r) • a = -v(r) e(r) • cr . (53) 
The TF-approximated Wigner function is then 

u(r, p) = V (-C - £g - V(r) + v(r) e(r) • <r) 
l + e(r) cr 



-v(p+{*)-p) 



1 — e(r) • cr 



r?(P_(r) -p) , 



(54) 



with 



P±(r) = [2M(-C-V(r)±v(r))} 1/2 , (55) 



and — C is the chemical potential. The underscore is a 
reminder that this Wigner function is 2 x 2-matrix val- 
ued. As a result, the single-particle density also has a 
corresponding spin dependence, 



l+e(r)-<T fP + (r) 
n(r) = tt' 



2nh 



l-e(r) a /P_(r) 

-7T 



2 " V ^nh 
|(n(r) + s(r) e(r) ■ cr) . 



(56) 



We observe that now there are two functions present here, 
the total density, n(r), and the spin- imbalance density, 
s(r), which are constrained by 



|s(r)|<n(r), 



(57) 



but are otherwise independent of each other. Therefore, 
the minimization to achieve the ground-state energy has 
to be done over both functions under the constraints of 
normalization and positivity: Eqs. (fT9|) and (|57| . respec- 
tively. 

We can then evaluate the trap, kinetic, and magnetic 
energy accordingly, 

Strap = tr 2x2 /(dr) \Mu 2 r 2 n(r) 



(dr) \Mu?r 2 n(r) , 



S km = tr 2x2 j {dr )-^ ^ u(r,p) 

= /(*)^[n(r) 2 + S (r) 2 ],, 
S m ag = -tr 2X 2 J{dr) n ■ B(r) n(r) 

'(dr) u(r)s(r) . (58) 



To compute the dipole interaction energy, we construct 
an approximation to the two-body spin-density matrix 
n( 2 )(r', r"; r', r") in the spirit of Eq. flSJ), starting with 
the single-particle orbital, 



^ x > ~ l J 8 m (x) 



(59) 



where a m and /3 m are the spin-up and spin-down com- 
ponents, x is the position variable, while x denotes the 
combination of the position and spin variables, so that 
the ground-state wave function of a N fermion system 
can be constructed as a Slater determinant, 



,x N ) = 



1 



det [(j> m {xi)] 
TV! m ^ 



(60) 
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When expressing the one-body and two-body spin-density matrices in terms of single-particle orbitals, we get 
n(x; y) = N / dx 2 ■ ■ ■ dx N ip(x, x 2 ,--- , x N )ip(y, x 2 ,--- , x N )* 



E 



a m (x)a m (y)* a m (x)/3 m (y)* \ = / n tt (x;y) n n (x;y) 



/3 m (x)a m (y)* /3 m (x)^ m (y)* 



n^(x;y) n u {x;y) J ' 



n (2) (x,y;x',y' 



N(N — 1) 



dx 3 ■ ■ ■ dx N ip(x, y, x 3 , ■ ■ ■ , x N )tp(x', y', x 3 , ■ ■ ■ , x N )* 



I . i n 



CCi(x) 

A(x) 



a m (y) 

My) 



az(x') 



am(y') 
/3 m (y') 



a™(x') 
/3 m (x') 



az(y') 



t(x; x') <g> n(y; y') - [n(x; y') ® n(y; x')] r 



(61) 



where T23 means interchanging the second and the third columns. The double summation in includes the I = m 
summands of the self-energy, which has equal contributions to the direct and exchange terms and hence does not 
contribute to the sum. 

The contact term in the dipole interaction potential is nonvanishing only in the singlet state, 



£dd, s = tr 4x4 J(dr')(dr") 



8tt 



S(r' 



L 



4tt 



(dr) [n(r) 2 - s(r) 2 ] , 



(62) 



where the prefactor y/2n' jl z originates in the reduction of dimensionality. We observe that, owing to the l/l z scaling, 
the relative strength of this term can be tuned by adjusting the stiffness of the z-confining trap. 

On the other hand, the triplet state interacts according to the remaining terms in the dipole potential. Since the 
state is symmetric under particle exchange, we use, instead of the in Eq. (|6T]) . an alternative two-body density, 



1 



n (2) (r', r"; r', r") = - ( n(r') ® n(r") - n(r'; r") <g> n(r"; r') 



(63) 



which yields the same energy but greatly simplifies the computation due to its tensor product structure. 
The triplet interaction energy is then given by 



-Edd.t — tr 4x 4 



47T 



{dv'){dv") n (2) (r',r";r',r") 



p 2 (7 ■ T — 3<T ■ pp ■ T 

? 



(64) 



with p = r' — r", and t denotes the Pauli vector for the second atom. To evaluate this expression, we apply the same 
procedure as that used to obtain Eq. ([29)) and find 



77 1 1 ' 32 ,/^T 1 / i 0/ i 



E. 



(i) 

dd,t 

(2) 
dd,t 



2 W 6 ^ 2 - 1 



(n(r) + s(r)) 5/2 + (n(r) - s{r)) 5/2 ~ ^(n(r) + s(r)) 3/2 (n(r) - s(r)) 



I /'0M 2 f {dr){dr ^ Vs z (r)-V >s z (r>) - V • s(r)V ■ s(r') 



2 4tt 



r - r' 



(65) 



where 7 = [P_(r)/P + (r)] 2 is essentially the ratio between the Fermi energies of the minority and majority spin 
components, and 



f(l) = (7- 1 + 14 + 7)^(7) + - 6 + 77)^(7) 



(66) 



is a combination of elliptic integrals that is smooth and finite for < 7 < 1. It is clear from Fig. [3] that 7(7) can be 
replaced by a linear function f(j) = ^j-tt + (16 — ^7r)7 to simplify computations. 

In passing, we note that a magnetic field with a large component in the xy-plane, such that 3e 2 (r) 2 — 1 < in Eq. 
(1651) for some region, may lead to an energy that is not bounded from below. The system is then instable and will 
collapse and explode within miliseconds (this catastrophe was observed in dipolar bose gases [13, [HI). When this 
happens, so much energy is made available that the 2D confinement will be lost. 
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For the simple case of a constant external magnetic field, we have 

B(r) = B Q e z , v(r) = B fi = v . 
In the dimensionless quantities, we re-parameterize the spin-imbalance density in accordance with 



(67) 



M x ) = 77 s ( r ) = cos (^( x )) ff( x ) 



(68) 



so that Eq. (|57p is automatically fulfilled. In the limit of weak interaction, we obtain the total energy as a functional 
of both g(x) and #(x), 



4fdM ] _ T 

HujN 3 / 2 



~ |(dx) 5 (x) 2 (l + cos 2 (tf(x))) +1 |(dx)x 2 5 (x) 



1 vq 

J^ 1 hcu 



(dx)0(x)cos(0(x)), (69) 



and the variation of g(x) and $(x), with \X 2 as the Lagrange multiplier, yields 



7rg(x) + cos(i?(x)) ( 7Tg(x) cos(i?(x)) 
g(x) sia(i?(x)) I 7rg(x) cos($(x)) 



I t; 



1 t>o 
TV 1 ftw 



0. 



(70) 



16 



/(7) 



1571-/4 



7 



FIG. 3: /(7) at relevant values of 7. It can be shown that 
15tt/4 < /( 7 ) < 16 for < 7 < 1. 



There are two nontrivial solutions, one for a spin- 
polarized (SP) system and the other allowing a spin- 
mixture (SM). The SP solution yields the TF profile 
equivalent to Eq. (|4o) . The SM solution gives 



with 



(2^/l ~ A 2 ' - x 2 ) for < 2: < z_ , 



1 

2^ 

— ( A + \J\ - A 2 ' - — ^ for < a; < x+ , 
I 27r \ 2 y 

(71) 



A = 



o/t 



(72) 



together with a constant spin-imbalance density in the 
center, 



,g(x) cos(i?(x)) 



A 



(73) 




FIG. 4: The dimensionless density profile of a SM (thin solid 
line) with A = 0.25, in comparison with the TF profile (i.e. SP 
solution, thick dashed line). The thin dashed line indicates 
the density of the majority component in the spin mixture. 
For greater values of A while keeping A < l/v 7 ?, the density 
profile approaches that of the SP solution with a lowering 
central density and an increasing x±. In the opposite limit, 
we recover a mixture of equal spin-components when there is 
no external magnetic field, i.e. A — 0. 



where x± are the radii of the spin-mixture (lower sign) 
and the entire cloud (upper sign) respectively, given by 



x± = 2(y/l-A 2 ' ±A) 



(74) 



These matters are illustrated in Fig. 2) 

However, the existence of a spin-mixture requires an 
extremely weak external magnetic field, such that 



1 1 



(75) 



This condition arises from the positivity of the radii 
x±. In usual experimental set-ups, this translates into 
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B ~ 10 mG for a system of N = 10 6 atoms with a ra- 
dial harmonic confinement of u> = 2ir x 20 Hz. In other 
words, a spin-polarized cloud is readily attainable. This 
justifies the treatment of a spin-polarized system before 
Sec. En 

VIII. SUMMARY AND OUTLOOK 

Table U summarizes the kinetic and the dipole-dipole 
interaction energies as functionals of the single-particle 
density for a fully spin-polarized gas in one, two, and 
three dimensions. It is clear that the structure of the 
density functionals depends crucially on the spatial di- 
mension. The procedure used here to reduce dimension- 
ality is by no means unique, but fairly well justified by 
the strong confinement of a stiff harmonic trap in a pos- 
sible experimental set-up. 

In 2D, the total energy with TFD approximation de- 
pends on both the dimensionless interaction strength e, 
and the number of particles N, as it does in 3D, but 
the A^-dependence is slightly more complicated. Namely, 
as one piece of the interaction energy is proportional to 
eA 1 / 4 , while the other piece is proportional to eA -1 / 4 , 
the latter is always a factor of ^/W smaller, inviting a 
perturbative treatment. 

For large eN 1 / 4 , the potential energy and the first piece 
of the interaction energy dominate, and are of the order 
of fiwiV 3 ' 2 x (eN 1 / 4 ) 2 / 5 , while the kinetic energy is of 
the order of HujN 3 / 2 x (eA^ 1 / 4 ) -2 / 5 . Numerical results 
suggest that the second piece of the interaction energy is 
of the order tkuN 3 / 2 x (eA^ 1 / 4 ) -158 , which makes it the 
most slowly growing term in the total energy. 

In the limiting case of e — > + , the kinetic and poten- 
tial energies are both fiwJV 3,/2 -y, where the equality is 



well predicted by the virial theorem applied to a simple 
harmonic oscillator. The two pieces of the interaction 
energy are of the order N 7 / 4 and N 5 / 4 respectively, even 
though both are vanishing due to small e. 

In addition to the fully spin-polarized situation, we also 
dealt with partially polarized gases, allowing for inho- 
mogeneous polarization. By considering the spin-density 
matrix, we found the energy as a functional of the to- 
tal density and the spin-imbalance density, and then de- 
termined the implied ground-state density profile. For 
typical experimental parameters and a modest number 
of trapped atoms, a spin-mixture can only exist for an 
extremely weak external magnetic field. In other words, 
a fully spin-polarized gas is readily attainable. 

Having thus established the TFD functionals, we in- 
tend to investigate the excitation energies of the sys- 
tem for small deviation from the equilibrium. On the 
other hand, it is well-known that the TF approximation 
is problematic at the boundary of the system. We will 
follow up on the gradient corrections of von Weizsacker 
type. It is perceivable that once the corrections are in- 

(2) 

eluded, may no longer be negligible. Lastly, we 

would like to explore other external trapping potentials, 
such as anisotropic harmonic traps, possibly with an op- 
tical lattice superimposed. 
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Appendix: Calculating the interaction energy 

The splitting of the dipole-dipole interaction energy into direct and exchange energy at the level of the one-particle 
density and the one-particle density matrix turns out to be inconvenient at it stands, since this yields two integrals, 
which both diverge individually but together sum up to a finite value. This prompts us to express everything in terms 
of the one-particle Wigner function, 



Eaa — — 



l^ 2 f(dr')(dr")(d Pl )(dp 2 ) 1 



2 4tt J {2nh) 4 |r'-r"| 3 

2 



Kr',piMr",P2) - ^(^, P iW^,p 2 )e 



3 i(pi-p 2 )-(r'-r")/fi 



(Al) 



r 



where we have used the substitution r=^(r' + r"), The integration over p can be evaluated with the outcome 
p = r' — r" and the Fourier transform in 2D, r i / 

/<*)£(<*»■ 

ik-i 



i-k 2 )-p/2 _ e i(p 1 -p 2 )-pfh 



Kr, p) = J (dk) „(k, p)e* - , = 2?r(|pi _ p2l/h _ |ki _ M/2) (A3) 

77(k, P ) = jM_ v (r,p)e-^ r . (A2) 
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TABLE I: Summary of the density functionals for the kinetic energy and the dipole-dipole interaction energy in one, two, 
and three dimensions. In ID, the spins are polarized normal to the z-axis, along which the atoms align. Erfc( ) denotes 
the complementary error function, and t = \z — z'\/(yTl±). Note that further simplification of the ID expression of E^d 
involves taking the limit of l± — > + which should be done with extreme care. In 2D and 3D, the spins are polarized along the 
z-direction, 6 in 3D refers to the azimuthal angle of the vector r — r'. 



-Ekin 



ID 
2D 
3D 




I dz dz ' ( n ( 2 ) n ( 2 ') - n ( z > z ') n ( z '\ z )j 



^ +t 2 )y^e z ~Ertc(t) -t 



(dr-0 



~^V^ n(ri) 5 ' 2 — 7rn(r)\/— V 2 n(rx) 



(dr)(dr) i n(r) - - ra(r') 



1 — 3 cos 



We recognize that the interaction energy is split into two 
pieces, 

^dd - ^ dd + ^ dd , 



E 



(1) _ I lion 



del 



M0M 2 o [,, \ (4>l) (jgg) , , W X 



2 (k 1+ k 2 ).r- (kijpi) - (k2jp2; 



iPl - P2 



(2) _ lMoM S 



dd 



2 4tt 



O /m N (^Pl) ( d P2) /JU w n \ 



Xe .( kl+k2 ). rF(ki)Pi) - (k2)P2) ^ik2| ) (A4) 



but this is not the splitting into the direct and exchange 
terms, as the integration of a single exponential term in 
Eq. (| A3|) will not converge. 

A closer look at E^ tells us that the J (dkx ) and 
J (dk.2) integrations recover the Wigner functions, which 
impose an upper limit of P = Uy/ 47rn(r) ' on the length 
of pi and P2 , so that 



E 



(i) _ 1 nop 



f^Ol^O f/A \ ^P 1 ) ( d P 2 ) I I /* 

— 2 ^ (dr) (2w(w lpl ~ p2l/ 



,M 2 4tt J ^''(2ttK) 2 (27Tft) 2 



4tt 



5/2 



(A5) 



(2) 

On the other hand, i? dd needs to be treated differently. 
The J (dpi) and J (dp 2 ) integrations yield the form factor 



e- lk r n(r) , (A6) 



while the J (dr) integration gives rise to a 2D Dirac delta 
function, 

J (dr) e i(kl+k2) - r = (2^) 2 5(kr + k 2 ) . (A7) 

This then takes care of one of the integrations over \s.i or 
k 2 , and we arrive at 



E (2) = _ l/fo^^a U dk )n( k )kn(-k) 



"dd 



2 4tt 



4tt 



tt /(dr) n(r)V-V 2 n(r), (A8) 



where the short-hand notation of %/— V 2 ' for the integral 
operator (|30|) is used in recognition that it is equivalent 
to the integral operator —V 2 when applied twice. This 
completes the derivation of Eqs. (f29|) . 
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